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Scroll waves are three-dimensional analogs of spiral waves. The linear stability spectrum of 
untwisted and twisted scroll waves is computed for a two-variable reaction-diffusion model of an 
excitable medium. Different bands of modes are seen to be unstable in different regions of param- 
eter space. The corresponding bifurcations and bifurcated states are characterized by performing 
direct numerical simulations. In addition, computations of the adjoint linear stability operator 

I eigenmodes are also performed and serve to obtain a number of matrix elements characterizing the 

■ long-wavelength deformations of scroll waves. 
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'. I. INTRODUCTION 

Depolarization waves in cardiac muscle, oxidation waves in the Belousov-Zhabotinsky (BZ) chemical medium or 
at the surface of certain metal catalyst and cAMP waves in colonies of slime molds are different examples of wave 
propagation in excitable media [0. They can be described in similar mathematical terms although the underlying 
\^ ■ processes are of a very different nature. In a two-dimensional (2D) or quasi-two dimensional situation, the propagation 
, of spiral waves has been observed in these three cases as well as in other excitable media Asides from their 

' intrinsic scientific interest, the potential role of these remarkable waves in cardiac arrhythmias and fibrillation |^ has 
motivated detailed studies of their properties during the past two decades. In particular, the mechanisms of different 
instabilities have been intensively investigated as well as their locations determined in parameter space of simple 
models and of experiments |^ . 

The potential relevance to cardiac dysfunction of scroll waves, the three dimensional analogs of spirals, has also 
been emphasized but their dynamics is still less thoroughly analyzed. Visualization of the BZ reaction in 3D gels 
^ I 1^,^ has confirmed the existence of these complex waves. Numerical simulations have revealed that they are prone 
to instabilities in several parameter regimes JlO|-p^. In order to analyze more systematically the different possible 
instabilities, we report here the result of computations of the full linear stability spectrum of a straight scroll wave 
O in a simple two-variable model of an excitable medium. This enables us to follow the different modes of deformation 
^ of a scroll wave and to investigate which type of modes become unstable in different regions of parameter space. 
^ The modes of the adjoint operator are also determined in order to compute the value of several coefficients given by 
matrix elements and to check proposed analytic relations. In addition, direct numerical simulations are performed to 
investigate the nonlinear fate of the different instabilities and to provide a detailed characterization of the restabilized 
' bifurcated states (when they exist). 

In section we define the studied two- variable reaction model and explain o ur n umerical methods. Some general 



> 



properties of spiral waves are also recalled. Then, we first consider in section HI a straight untwisted scroll wave. 
It is the simplest (i.e. z-independent) extension in the third (z-) dimension of a two dimensional (xy) spiral wave. 
At the linear level, two possible types of instabilities are found. Modes with positive real parts can be observed on 
the translation bands, which correspond to z-dependent translations of the 2d spiral in the different xy-planes, or on 
the meander bands, which come from z-deformations of the 2d spiral meander modes. At the non linear level, the 
translation band instability gives rise to a scroll wave with a continuously extending core and does not lead to 

a restabilized non linear state. We confirm that this type of instability is directly related||l4| to the drift direction 
of a 2d spiral in an external field. In contrast, the meander band type of instability JlSf generally restabilizes in a 
distorted scroll wave and no simple relation to 2d spiral drift is observed. 



In section IV, we consider twisted scroll waves. A 3D steady wave is built by rotating (i.e. twisting) the 2d 



spiral around its rotation center as one translates it along the z-direction. We find that twist exceeding a definite 
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FIG. 1. Spiral bifurcation diagram for Eq. with e — 0.025 and the reaction terms of ref. [19]. The bold line dM is the 
meander threshold instability line and separates steadily rotating spirals (SS) from meandering spirals (MS) (above b = 0.02 the 
thin line denotes our less accurate determination of dM). The line dR marks the boundary of spiral wave existence on the left 
of OR the wave tip retracts (see e.g. ref. [14] ). The symbols along the line b — 0.01 denotes the different parameters at which 
scroll waves are studied in the present work. Stars represent stable scroll waves, crosses (x) and pluses(+) meander-unstable 
scroll waves, circles (o), and (•) translation-band- unstable scroll waves. 2D meandering spirals are represented by (-I-) and (o). 



threshold can lead to the appearance of unstable modes in the translation brands of the scroll wave. This "sproing" 
pof instability is seen to take place a finite wavevector away from the scroll wave translation symmetry mode. We 
provide analytical arguments which show that this very generally results from 3D rotational invariance in an isotropic 
medium. Nonlinear development of this instability when a single unstable mode is present results in a restabilized 
helical wave, as previously described [Q. Properties of these nonlinear states are computed and compared with the 
linear characteristics (wavevector, frequency) of the sproing instability. When several unstable modes are present in 
the simulation box, the scroll wave core filament takes a more complex shape which is found to travel like a nonlinear 
wave of constant shape in the vertical direction. Three appendices provide the details of our numerical algorithms, a 
derivation of a general formula for spiral drift in an external field and useful ribbon geometry formulas. A fourth one 
explains the relation between the present calculations and previously derived averaged equations for the motion of a 
weakly curved and twisted scroll wave . 

The results of our linear stability analysis have previously been shortly described in p7|]. 



II. METHODS AND GENERAL RESULTS 



A. Reaction-diffusion model 



Two- variable reaction diffusion systems have been shown to describe semi-quantitatively spiral waves dynamics and 
its "generic" features. They have been used in various contexts since their original introduction as a simplification 
of Hodgkin-Huxley dynamics. The analysis of such a simple model appears in any case as a useful first step before 
going to a more complicated description if required. We thus follow this classic path and take for the excitable medium 
dynamics 

dtU = V^u + fiu,v)/e (1) 
dtv = g{u,v) (2) 

We only consider the singly diffusive case, the most relevant to cardiac physiology. For definiteness, we also choose 
reaction terms /(m, u) = w(l — u)[m— (u-|-6)/a], g{u^v) = u — v as proposed in [p9[. This permits fast direct simulations 
and tests of our numerics by comparison with previous results for spiral waves. The 2d-spiral bifurcation diagram 
for this model is shown in Fig. |^ for variable values of the parameters a and & at a fixed value of e = 0.025. The 
meander instability line {dM with the notation of [^) is plotted with its "large core" branch at smaller values of a 
than its "small core" branch. The crossing line which separates meander trajectories with outward petals from those 
with inward petals is also drawn as well as the diverging core existence boundary of spiral waves (OR). 

In the following, the stability and dynamics of scroll waves are analyzed at several points along the line b = 0.01, 
as a varies and crosses the different boundaries. 
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B. Numerical strategy 



In order to obtain the steady scroll waves and compute their stability spectrum, Eq. (|l|, ^ are written in cylindrical 
coordinates with u and v functions of r, (j) = 9 — cot — z , t and z, 

[dt + 2T^dl, - dl,)u = {ud^ + rldl^ + VIj,)u + /(u, v)/e (3) 
dtv = ujd^v + g{u, v) (4) 



1. Steady states 

A steady scroll wave for a given imposed twist Tw is a time independent solution of Eq. ^) with u{r, (f>, t, z) = 
uo{r, (j)) and w(r, 0, t, z) = vo{r, 0) and rotation frequency w = cui, 

{VId + ^id<f, + Tldl^)uo + /(wo, vo)le - (5) 
uiid^vo + g{uo,VQ) = Q (6) 

This nonlinear fixed point problem for the function uq^vq and n onlin ear eigenvalue oji is solved by using Newton's 
method after discretization of Eq. (^,^, as detailed in appendix Al. It should be noted that Eq. (||, ^) are purely 



two-dimensional due to the scroll wave translation symmetry along the z axis in the introduced coordinates. 



2. Linear stability 

Once a steady scroll wave is obtained, one can linearize Eq. ^ around it. Invariance of Eq. ^ ^) and of the 
steady state by translation along the z-direction (in the introduced twisted rotating frame) leads to decompose a 
general perturbation depending on the three spatial coordinates on its Fourier components along the z-axis. We thus 
consider perturbations under the form u — uq + eKp[a{kz)t — ikzz]ui{r, (j)),v = vq + exp[a{kz)t — ikzz]vi{r, cj)). The 
linear equations obeyed by and the (complex) growth rate a{kz) read, 

aui = {-kl + 2iTy,kzd^)ui + (ujid^ + r^c)^^ + ^Id)'^i + [5«/(wo, wo)wi + 5i,/(uo, wo)wi]/e (7) 
avi = wi^^ui + dudiuQ, vo)ui + dvdiuo, vo)vi (8) 



or symbolically. 



So in a second step, the (~ 10) eigenvalues of largest rea l pa rts of £k^ are precisely determined (for given kz and 
Tji,) using an iterative algorithm [ pO[ detailed in appendix |A 2| . The whole numerical procedure is quite similar to 
the spiral linear stability analysis of rcf. and extends it to 3D. It should be noted that the twist rate can be 
prescribed at will so the procedure is not confined to weak twist (of course, for too large a twist, a steady scroll may 
no longer exist and the Newton steady state finding procedure fails to converge). 

Two points are worth emphasizing: 
-taking a Fourier transform has eliminated the z-direction so Eq. (0, |^) are purely two-dimensional (but fc^-dependent) 
as the steady state equations (||, ||), 

-correlatively, each mode of the 2D-spiral is replaced by a band of modes indexed by the wavevector kz- 
Some general properties of the spectrum can be noted at this stage. 

For zero twist, Eq. (0, ||) depend only on k^ so the spectrum bands are even functions of kz- Moreover, Ck^ is a 
real operator and its complex eigenvalues come in complex conjugate (c.c.) pairs. So, bands of complex modes also 
come in complex conjugate pairs. 

For non-zero twist {tw ^ 0), these symmetries non-longer hold. It only remains true that Ck^ — C*~k^- So bands 
of complex modes can be grouped in pairs <Ji{kz),(J2(kz) with 172(^2) = <yi(—kz) 
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3. Direct numerical simulations and instantaneous filament definition 



In order to determine the bifurcation type a nd th e fate of each instabihties, we performed direct numerical sim- 
ulations of Eq. ([^,|^) as explained in appendix A3. In two dimensions, it is usual to define the spiral tip as the 



intersection point of two (somewhat arbitrary) particular level lines u = Utip,v = vup. The spiral tip trajectory is 
then a convenient way to visualize the spiral dynamics and its core instabilities. Similarly in 3D, we choose here 
to define the instantaneous filament as the intersection line of the two particular level surfaces u = uup = 0.5 and 
V = Vtip = 0.75(0. 5a — b). It can be thought of as the line of spiral tips in the different xy-plane. 



C. Special eigenmodes 

In the 2D case, there are five dominant modes of spiral dynamics in the simple case described by Eq^ (0, |) [|2) : 
-one neutral mode with cr = 0, the rotation mode, which comes from the rotational invariance of Eq. (y, 
-two complex conjugate purely imaginary modes with a = iicji, the translation modes, coming from the translation 
invariance of the starting equations (|l|, |^) 

- two complex conjugate modes, the meander modes, corresponding to the oscillatory meander instability the real 
part of which crosses zero on the meander instability line. 

In the following, the five bands of modes originating from these special modes are found to play the most important 
role in scroll wave dynamics in the sense that they have the largest real parts and that each instability of a scroll wave 
can be ascribed to one of these bands (i.e. the modes on a part of that particular band acquire a positive real part). 

Thus, before proceeding it is worth recalling the expression of these symmetry eigenmodes for spirals pl| and their 
straightforward generalization for scroll waves. 

The rotation mode is the simplest. Differentiation of Eq. (|^, ^) directly shows that (9^mo, cJ^fo); the rotation mode, 
is a solution of Eq. (0, ||) for = Q and tr = 0. 

The eigenmodes associated to invariance under x, y-translations of the spiral rotation center are less straightfor- 
wardly obtained because the spiral is a time independent solution in a rotating frame (with (x',y') coordinates). A 
steady spiral rotating around the origin is given by 

Uo[x' , y'] — Uo[cos{u}it)x + sin{ujit)y, — sm{ujit)x + cos{uJit)y] (10) 

The corresponding spiral rotating around the point (xq, j/q) is then 

Ua[cos{ujit){x - xo) + sm{ujit){y - yo), - sin{ujit){x - xo) + cos{uJit){y - yo)]. 

Using the {x',y') rotating coordinates, this translated spiral reads, 

Uo[x' — Xq cos{uJit) — yo sin(cjit), y' + xq sin(ci;ii) — j/q cos(cL'ii)] = 

Uo[x -( ^ exp(iwit) -|-c.c.),y ~ (i ^ exp(«wit) -f c.c.)J (11) 

Expansion of Eq. ( pT| ) for small XQ,yQ gives the expression of the translation modes 

f^A^f i f + ^dy'l^o \ ^ ( d uo + ^^,uo/r \ 

\Vt J \ [Ox' +tdy')Vo J t-y y drVo + ld4,Vo / r J 

with eigenvalue iwi, and the complex conjugate eigenvector with eigenvalue —iuji. 

This computation immediately generalizes to untwisted scroll waves. This is also the case for twisted scroll waves 
but one should recall that Eq. (|^, |^) are written in a referential that rotates in time but also as one moves along the 
z-axis. This modifies the exponential factor in Eq. ([ll| ) which becomes exp(ia;ii -I- t^z) to include the z-rotation. So, 
for a twisted scroll wave the translation mode (ut,Vt) is an eigenvector of Ck^=-Ti^ with eigenvalue ilui. The other 
(complex conjugate) translation eigenvector (u^^^v^) is associated to the eigenvalue — zwi of the now different linear 
operator Cu^=r^ . 

A direct algebraic proof of these facts can be given. If ones defines the two operators 

T = exp(i0)(9r + i/rd^) 

M = ujid^ + rldl^ + V^^ (13) 
a direct computation gives the commutators 
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[T, M] = -iuiT - r2 (1 + 2id4,)T 

[T,d^]=-iT (14) 

With these notations, the fixed point equations reads Muq + /(uo,wo) = 0, LUid^vo + g{uo,vo) = 0. Action of T on 
these two equations gives TMuq + duf Tuq + dvfTvo = 0, uJiTd^VQ + dugTug + d^gTv^ — 0. Commutations of 
T and M (in the first one) and T and 9^ (in the second one) using Eq. |l^ directly show that {ut,vt) = {Tuq,Tvq) 
satisfies Eq. (^ ||) with cr = +iuii and = —t^,. 

D. Left eigenvectors and scalar product 

The Unear stabiHty computation can be extended to determine the left eigenvectors of Ck, We have found it worth in 
particular to compute the left eigenvectors of the 2D spiral stability operator C (i.e. Ck^=o for = 0) corresponding 
to the translation and rotation modes since they often appear in perturbation calculation [ |l6yi2|] (for examples, see 
III Al and appendix^). 



section 



We simply define the scalar product between a left {u,v) and right {ur,Vr) two-component function by integration 
over the whole 2D space as {u, Ur) + {v, Vr) where 



{f,g) = JJdrd^rfir,^) gir,^) (15) 

The left eigenmodes of C thus obey, 

aui = {-uJid^ + V2£,)mi +duf{uo,vo)ui/e + dugiuo,vo)vi (16) 
avi ^ -uJid^vi + dyf{uo,vo)ui/e + dyg{uo,vo)vi (17) 

The result of one such computation for the left translation eigenmode {utjVt) (the solution {ui,vi) for a = iuji) is 
shown in Fig. |^. In contrast to the translation mode {ut,Vt), the left eigenmode {ut,Vt) quickly decreases away from 
the spiral core as argued in ]l^ and explicitly obtained in the free-boundary limit ||lj,^ (but opposite to what is 
supposed in [0). This also holds for the left rotation mode {u^,v^) (the solution {ui,vi) for a = 0) as shown in 
Fig. ||. As a consequence, the scalar product (^ between these left functions (m, v) and any right function (u^, v^) 
(even slowly increasing) is well defined^. We do not find it useful to include a time integration in the scalar product 
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III. UNTWISTED FILAMENTS 

We begin with the simplest case of scroll waves with no twist (r^, = 0). In this case, the steady scroll equations 
(^, P) are clearly identical to those of a 2d spiral. We take as an example the parameter value a — 0.9, b — 0.01 (and 
e ~ 0.025). A steady spiral/scroll wave is found for a rotation frequency loi = 1.769. The linear spectrum of modes 
around this steady scroll is plotted in Fig. ^ (only the upper quadrant upper fc^ > 0, Im{a) > is shown since the 
other quadrants can be deduced by parity and complex conjugation) The five translation, rotation and meander modes 
of the spiral wave stand at fc^ = 0. The steady spiral is stable as shown by the negative real parts of the meander 
modes. As stated above, the spectrum around the scroll wave is organized in several bands of modes which originates 
from the spiral modes at kz = 0. Only the five less stable bands are shown in Fig. ^. At these parameter values, 
extension to the third dimension does not bring any instability (at least at the linear level) since as seen on Fig. ^ 
the real part of cr(fcz) becomes more negative on each band as increases. For other parameter values, a straight 
scroll wave can however be unstable while 2D-spiral are stable. This can happen in two different ways. Depending on 
position in parameter space, either the translation or the meander bands become unstable for k^ ^ 0. We examine 
these two cases in turn in the following two subsections. 



^The fast decay of the left eigenmodes makes the space integration converge without any additional factor. Adding one such 
extra factor as suggested in pii| would actually make all the scalar product vanish. 
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FIG. 2. Contour plots of the translation eigenmode, (a) u-component modulus, (b) v- component modulus, and of the 
corresponding left eigenvector , (c) u-component modulus, (d) v- component modulus. The maximum value of the fields is 
set independently for u and v equal to 1. and the contours are plotted for {a.)u = 0.01, 0.05, 0.1, 0.3, 0.4, 0.5, 0.6 and 0.7, 
{h)v = 0.01, 0.05, 0.1, 0.2, 0.4, 0.6 and 0.8, (c) and (d) u,v = 0.0001, 0.001, 0.01, 0.05, 0.1, 0.3, 0.5, 0.7 and 0.9. The 
parameter values are a — 0.44, b = 0.01 and e = 0.025. The pulsation of the steady rotating spiral is a; = 1.1612. The circles 
represent the limit of the simulation box. 



A. Translation band instability 



The translation bands can have unstable modes for while the 2D spirals are stable. An example of this 

phenomenon is shown on Fig. ^ for a = 0.44, b = 0.01 and e — 0.025. A qualitatively similar spectrum is obtained for 
all points in Fig. |l| denoted by filled dots (•) near the large core spiral existence boundary dR. As seen on Fig. |, the 
instability takes place for small as soon as kz is non-zero. It corresponds to the "negative line tension" instability of 
ref. . We show below that the curvature of the translation modes at fc^ = is given by the spiral drift coefficients 
in an external field. So, this translation band instability is directly related to the fact that 2D spiral drift opposite to 
an applied external field in this parameter region. 



1 . Long wavelength stability and 2d spiral drift in an external field 



As recalled and shown in detail in appendix a small applied external field E induces a drift of the spiral rotation 
center at a velocity v such that 

V = a\\E + a^iiJi X E (18) 

where uji is the spiral rotation vector. It has previously been noted that a weak scroll wave curvature acts as an 
external field and therefore that a straight scroll wave is unstable if a|| < since a small curvature tends to grow. 
More precisely, the small kz behavior of the two translation bands is given by 

cr±{kz) = ±iuJi + ± ia±) kj + 0{kt) (19) 

Eq. (19) is simply derived by a first-order perturbative calculation as follows. The linear eigenvalue problem (^, |^) 
reads 
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FIG. 3. Contour plots of the rotation eigenmode, (a) it-component modulus, (b) v- component modulus, and of the 
corresponding left eigenvector , (c) u-component modulus, (d) v- component modulus. The maximum value of the fields is set 
equal to 1. and the contours are plotted for (a) u = 0.001, 0.1, 0.3, 0.5 and 0. 7 , {h)v = 0.001, 0.1, 0.2, 0.4, 0.6 and 0.8, (c) 
and (d) u, v = 0., 0.001, 0.01, 0.1, 0.3, 0.5, 0.7 and 0.9. Same parameters as in Fig. g. 
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FIG. 4. (a) Real and (b) imaginary parts of the growth rate a{kz) as a function of the wavenumber fcz for the translation 
(+), rotation (•) and meander (o) bands. The parameter values are a = 0.9, b = 0.01 and e = 0.025. Fig. ^ shows that 
the meander mode at fez = is stable and that the growth rate decreases on the meander band with k^. The translation 
mode is also restabilized for finite values of fcz. The translation and meander bands are well approximated by respectively 
a-t(fcz) = il.769-|-(-0.65 + j0.61)fcf and cr„(fcz) = -0.3411 -|-il. 720 + (-0.25 -|-j0.01)fe^ The value of at (0) is in good agreement 
with the independently determined pulsation of the steady scroll wave oji — 1.769 
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FIG. 5. (a) Real and (b) imaginary part of the growth rate a{kz) as a function of the wavenumber for the parameter value 
a = 0.44, b = 0.01, e = 0.025. The meander mode, (o) is stable and the growth rate decreases with the wave number on the 
meander band. The modes of the rotation band (•), are also stable. The mode of the translation band (+), are unstable for 
finite values of fcz. 
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For small kz, the modes of the translation bands can be obtained by perturbation around the known translation 
modes at fc^ = (Eq. p^ ). For definiteness, we consider the upper band (which start at a{kz = 0) = +iuji) and write 
a{kz) — iLUi + 5a where (5cr ^ 1 is the sought perturbative correction, 
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Substitution in Eq. (E(H) gives 
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The first-order expression of 6a is obtained in a usual w ay by taking the scalar product of Eq. 
eigenvector {ut,vt) of Ck^=o for the eigenvalue icui (sectionilD) 



(21) 



(22) 



m) with the left 



da 



{Ut,Ut) 



{ut,ut) + {vt,vt) 



(23) 




since the matrix coefficients on its r.h.s also gives the spiral 



Eq. (|2^) is equivalent to the announced formula 
drift coefficients, as shown in appendix ^ (sec Eq. 

In Table ^, the spiral drift coefficients a|| and a± are compared to the results of independent computations of the 
curvature of a{kz) at fc^ = for the translation bands, from diagonalisations of Ck^ at different values of a. The good 
agreement between these results is a check both of the analytic formula ( |l9| ) and of our numerics. To recapitulate, 
the translation band instability is found to be a long wavelength instability (i.e. the band of unstable wavelength 
starts at kz = 0) which is present in the whole domain of parameters where a 2D-spiral drifts opposite (given our sign 
convention in Eq. Bl) to the applied field. 



2. Nonlinear evolution of the instability 

The nonlinear fate of this translation band instability was studied by direct dynamical simulations of Eq. (^, ^) . In 
the parameter regime of Fig. ^ when modes of the translation bands are unstable for finite values of kz, an initially 
straight scroll wave was observed to be unstable provided that the simulation box was large enough to accommodate 
an unstable mode. In agreement with previous observations |l2| , the filament was observed to increasingly depart 
from its straight initial configuration and its length was observed to grow in the simulation box. When the filament 
eventually collided with the boundaries of the simulation box, it split into two filaments. This repeated again and no 
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a 




a['{k, = 0)/2 


a^k. = 0)/2 


ay — ia± 


0.44 


1.16 


1.9 + 0.82i 


-1.6 + 0.78i 


-1.97 - 0.84j 


0.48 


1.38 


3.2 + 0.44i 


-3.7+ l.lOi 


-3.0 - 0.49j 


0.67 


1.76 


-2.14 + 0.85i 


1.61 + 0.25i 


2.2 - 0.9j 


0.7 


1.78 


-1.63 + 0.83i 


1.04 + 0.2H 


1.62 - 0.83i 


0.8 


1.81 


-0.87 + 0.701 


0.08 - 0.06i 


0.854 - 0.7H 


0.9 


1.77 


-0.65 + 0.61i 


-0.25 - 0.089i 


0.66-0.61i 



TABLE I. The scroll wave pulsation uji, half the second derivative a"{kz — 0)/2 of the translation band at kz — (with 
a"t(0) = iuJi), half the second derivative a-'^{kz — 0)/2 of the meander band at fc^ = (with lm{am{0)) > 0) and the drift 
coefficients of the 2D-spiral in an electric field Qfy — ia± for b = 0.01, e — 0.025 and different values of a. 






FIG. 6. Instantaneous filament evolution starting from a slightly perturbed straight scroll for equally spaced times 
{t — 25., t = 50., t = 75. and t — 100.) during a simulation in a simulation box of size (128 x 128 x 120) with a space 
discretization step dx = 0.2 using periodic boundary conditions. The parameters are a = 0.44, h — 0.01 and e = 0.025 and 
correspond to the linear spectrum shown in Fig. ^. 
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restabilisation was observed. A typical evolution is shown in Fig. ^. The minimum simulation box size which allowed 
the instability development closely agreed with the results of the linear stability analysis. For instance, in the case 
a = 0.44, b = 0.01, and e = 0.025, the maximal kz of the unstable band is obtained to be kz = 0.84 via the linear 
stability analysis whereas the direct numerical simulations show a minimal simulation box height corresponding to 
kz = 0.81 

The choice of boundary conditions on the simulation box top and bottom faces influences the minimum box height 
for the instability development (but, apart from that, was not observed to qualitatively modify the instability nonlinear 
development). This critical size was found to be twice bigger for periodic periodic boundary condition than for no- flux 
boundary conditions which can accommodate linear modes of wavelength twice as long as the box height. 

B. The third dimension-induced meander instability 

In the region where a 2D-spiral drifts toward an applied field (i.e. a|| > 0), the modes of an untwisted scroll wave 
translation bands are stable. As pointed in [^ , an untwisted scroll wave can nonetheless be unstable in a parameter 
region where a 2D-spiral is stable. This happens when the meander bands are destabilized by deformation in the 
z-direction as we study below. 

1. Linear analysis 

This induction of the meander instability by three-dimensional effects is shown in Fig. ^ For the parameters of 
Fig. I^a, all modes have negative real parts and the scroll wave is stable. However, one sees that the real part of 
the modes on the meander band starts by increasing as kz increases from zero. For the parameters of Fig. ^ which 
stand closer to the 2D meander boundary, a finite band of modes with kz ^ has acquired a positive real part while 
the real part of 2D spiral meander mode at = is still negative . Thus, for these parameter values close to the 
" small core" side of the 2D meander instability boundary, a 3D scroll wave is unstable to meander while the steadily 
rotating 2D spiral is still stable as pointed out in . On the " large core" side of the meander instability boundary. 
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FIG. 7. The growth rate Re[a{kz)] as a function of the wavenumber for the parameter value (a) a = 0.7 and (b) a — 0.67(b) 
The other parameter are b — 0.01 and e = 0.025. The meander mode, (o), is stable for fcz = both in cases (a) and (b). The 
growth rate increases on the meander band with kz . In case (a), it always remains negative and there is no instability. In case 
(b), it becomes positive for fez higher than 0.30 and lower than 0.69 showing the finite-fc^ instability of the steady scroll wave. 
The rotation (•) and translation bands (-I-) are stable. In the phase diagram Fig.0 the points with a spectrum similar to the 
(b) case are represented by crosses (x). 



three-dimensional modulations have on the contrary a stabilizing effect on the meander instability and the 2D and 
3D mea nder instability thresholds coincide as shown on Fig || (note however that the translation band instability of 
section [II A renders the scroll wave unstable for these parameters) . The fact that the curvature of the translation 



■^For the observed largest unstable k^, numerical simulations were not carried long enough to observe the full nonlinear 
development of the instability. However, it was observed with a box height increase of a single space step. 
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(b) 




FIG. 8. Real and imaginary part of the growtli rate a{kz) as a function of the wavenumber for b = 0.01, e = 0.025. (a), (b): 
a = 0.5; (c), (d): a = 0.52. In both cases, the real part decreases on the meander band and increases on the translation band 
with fez for small values of kz- In the a = 0.52 case, hybridization of the eigenvectors of the two bands exchange these trends 
at slightly higher values of fez i.e. the real part of the meander band increases whereas the real part of the translation band 
decreases. In the phase diagram of Fig. 0, the points where the linear spectrum is similar to one of these cases (with both the 
meander and translation bands unstable at small fez) are represented by represented by circles (o). 



and meander bands are of opposite sign in Fig. |5| and Fig. [7| may lead to think that the finite behavior of the 
meander bands is also related to 2D spiral drift (as indeed proposed in p^). However, a quantitative computation 
shows no simple relation between the meander bands curvature at — and the spiral drift coefficient as reported 
in table |[ Moreover, even at the qualitative level, there is no general validity to the opposite sign rule between the 
translation and meander bands curvature at fcz = 0, as shown by the data of Fig. ^. 

It is of some interest to see how the spectrum of Fig. |^ is transformed into the spectrum of Fig. |^ as one traverses 
the 2D meander unstable region. This happens through hybridization between the translation and meander bands as 
illustrated in Fig. ||. 



2. Restahilized bifurcated states 

We performed direct numerical simulations to examine the nonlinear development of the meander bands instability 
at fcz 7^ and to characterize the restabilized nonlinear states. For boxes of medium size in the z-direction as used 
in our numerical simulations, the restabilized nonlinear states strongly depend on the top and bottom boundary 
conditions. Two different ones were implemented (no flux boundary conditions were always implemented on the 
side boundaries). We used either periodic boundary conditions which permit the development of a single pair of 
unstable modes and are simpler to analyze or no-flux boundary conditions since they are clearly more relevant in an 
experimental context. We describe the results in turn. 

Periodic boundary conditions. 

Direct numerical simulations were performed in the parameter regime of Fig. |^. Three types of initial conditions 
were used which all led to a stationary restabilized scroll wave after a transient regime. 

For the first one, the perturbation of the steady state was chosen to contain the unstable wavelengths ifc^. The 
initial fields were chosen in the form u{x, y, z) = U2d{x, y)(l + a cos(fczz)), v{x,y, z) = V2d{x^ y){\ + a svaQtzz)) with 
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the magnitude a of the 3D perturbation being of order 10~^ (it2Z3 and V2d correspond to a steady spiral wave which 
is stable in 2D in this regime) . This perturbation transforms the straight unperturbed instantaneous filament into an 
"helix" of small elliptical cross section and pitch 2-K/kz. When the wavelength kz corresponded to an unstable mode, 
this helix was observed to grow and to reach a helical restabilized state of periodicity k^. 

A slightly more complicated time development was observed with an initial perturbation of the straight scroll wave 
in the form u{x,y,z) — U2D{x,y){l + acos{kzz)),v(x,y, z) = V2Dix,y){l + a cos(fczz)). This gives the instantaneous 
filament {x f (z) , y f (z)) a planar "zig-zag" shape of small amplitude, x{z) — xi cos{kzz), y{z) = yicos{kzz). The 
zig-zag perturbation was first observed to grow before turning into an helical restabilized state of periodicity kz ■ This 
type of time development is expected on general ground in system where left and right progressive wave compete (see 
e.g. 11). 

Finally, competition between several different unstable modes (up to 4 different kz) was examined with initial 
conditions such as u{x,y,z) = U2Dix,y)[l + aexp(-(z - Zo)"^ /Lc)],v{x,y, z) = V2D{x,y)[l + (3exp{-{z - ZoY /Lc)]. 
(typical values are ||(a,/3)||2 = 0.01 and Lc about a few tens of space steps). In that case, it was generally observed 
that the most unstable wavelength compatible with the box height was selected. At a qualitative level, the transient 
regime was found to mainly depend of the planar or non planar character of the initial perturbed filament : When it 
was planar (case (3 = a) & zig-zag filament first grew before taking an helical shape. When it was non planar, [3 ^ a 
then an helical filament grew directly. 

Close to the instability threshold, the shape of the restabilized instantaneous filament is closely approximated 
by an helix of circular cross section, as shown in Fig. |9[ and its motion can be portrayed in a way similar to the 
epicycle description of meander. The axis of the instantaneous filament helix rotates around a fixed vertical axis at 
the frequency of the steady two dimensional spiral. In the rotating frame where these two axes are fixed, the helical 
instantaneous filament itself rotates with a frequency close to the imaginary part of the meander linear eigenmode 
(i.e. the difference between these two frequencies is of order 10~^ in a strong meander regime and of order 10"'^ in a 
weak meander regime). Thus, for periodic boundary conditions the meander amplitude is found to be independent of 
height (z) while the phase of the epicycle motion varies linearly with z. 

We performed two different systematic studies in order to better characterize the nature of the 3D meander bi- 
furcation. In the first one, we kept the size of the simulation box constant and varied the excitability using a. In 
the second one, we kept a constant and varied the size of the simulation box, that is the wavenumbcr of the initial 
perturbation. 

At the linear level, the results of these direct numerical simulations are in close agreement with the predictions of 
the linear stability analysis, both for the the instability threshold ac (with an accuracy of order 10"'^) and for the 
unstable wavenumber range [fc_(a), k+{a)]. 

The radius of the instantaneous filament helix can be taken as a measure of the 3D meander instability amplitude. 
As shown in Fig. it is found to behave as the square root of the distance to the threshold (either |a — ad or 
\kz — k±{a)\). Therefore as reported previously p^ ], the 3D meander bifurcation is a supercritical Hopf bifurcation, 
as in 2D. 



However, the kz dependence of the nonlinear term quickly becomes important away from threshold. In Fig. 10, the 
square of the meander amplitude is compared to the growth rate of the unstable meander mode in the whole band 
[fc_(a), k+{a)]. A clear asymmetry of the curve is already seen, with a slope at the /c_(a)-end about 3.4 times 
larger than the /c+(a)-end slope. 

No flux boundary conditions. For the box height H used in our simulation (about 2 or 3 unstable wavelengths), the 
boundary conditions chosen on the top and bottom boundary conditions have a strong influence on the restabilized 
state. For no-flux boundary conditions, the spiral wave in each horizontal xy— plane has a meandering motion. 
However, contrary to the case of periodic boundary conditions the meander amplitude depends on z and is well 
approximated by |cos(fcz2;)| (see Fig. This implies in particular that in in some horizontal planes the meander 
amplitude is zero and that the corresponding spiral tip performs a simple steady rotation. In the rotating frame 
where these special tips are motionless, the instantaneous filament takes a planar zig-zag shape that rotates around 
its vertical midline at a pulsation close to the imaginary part of the unstable meander eigenmode corresponding to 
the wavelength of the filament. Thus, for no- flux boundary condition, the meander amplitude varies sinusoidally with 
height while the phase of the epicycle motion is independent of height. 

This shape and motion are simply understood in a weakly nonlinear description where the restabilized state can be 
approximated as a sum of the unperturbed solution and the four unstable meander eigenmodes {(Jm{kz) and a^{kz) 
at ztkz) 

u = uo{r,tP) + [Auiir,i}) e'('=^^+'^^*) Bui(r, V) e'^'''''''+'^''^ + c.c] (24) 
V = Vo{r,tP) + [Aviir,^;) e^('=-"+"^*) Bi;i(r, V) 6^^^'''^'+'^''^ + c.c] (25) 

with ui, vi the eigenmode of eigenvalue am{kz) and uj2 = Im[am{kz)\ the meander frequency. The no- flux boundary 
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condition , dzU = 0, at z = and z = H = 27r/fc^ enforces A = B (i.e. for the height considered, the no-fiux boundary 
conditions stabiUze the state with symmetric upward and downward propagating deformation which was observed 



to be unstable with periodic boundary conditions). The instantaneous filament corresponding to the fields (25) is 
easily determined by remembering that it is the locus u — Wtip, v = Utip. Its position is conveniently parameterized as 
a;jip(z) = a;Q + (5x'(z, t), y[-^^{z) = yQ + Sy'{z, t) using Cartesian coordinates in the rotating frame where the unperturbed 
filament is standing at (a;g,?/g). For small \A\^ one obtains 

d^oiQ 5x' + dyOiQ Sy' + cos(A:^z)[2Auie''^"* + c.c] = (26) 
d.j:a}a Sx' + dyOiQ Sy' + cos{kzz)[2Avie^'^^^ + c.c] = (27) 

where the field ui, wi, uq, vq and their derivatives are evaluated at the unperturbed filament position {x'q, y'^). Inversion 



of Eq. <m, 27) gives 



Sx' = cos(fc^z)[aAe^'^=* + c.c] 

Sy' = cos{k,z)[(3Ae''^^^ + c.c] (28) 

where a and [3 are complex constants which depend on ui, vi and the derivatives of mq evaluated at the point (xq, yo). 
This clearly shows the planar zig-zag shape of the instantaneous filament since points in different xy-planes simply 



differ by the real scale factor cos{kzz). As time evolves, Eq. |28| also shows that the filament points follow scaled 



ellipses in different xy-planes. Our simulations show that, as for 2D meander, these ellipses are in fact almost circular 
0- . . . 

Finally, we briefly discuss the transition between the restabilized meander regime seen on the "small core" side 
of the phase diagram (Fig. |l|) and the negative line tension dynamics which belongs to its meander "large core" 
side. Since scroll waves do meander in this transition region, the evolution of scroll waves as seen in direct numerical 
simulations is not directly linked to the linear spectrum of the steady scroll wave. For example, when the meander 
and translation bands are strongly hybridized, the translation bands are only unstable for small values of kz (Fig. |^). 
Nonetheless, direct simulations show that the meandering scroll is unstable and that its core grows as in the negative 
line tension regime in simulation boxes small enough to only contain unstable modes of the meander bands with larger 
values of kz ■ This happens on the whole small a side of the dashed uii = ui2 line of Fig. |l| This line stands very close 
to the line where the external field drift of meandering spiral changes sign and it is difficult to distinguish the 
two in our simulations. So, the link between the spiral drift sign and the "negative line tension" type of instability 
development continues to hold for meandering scroll wave. 



IV. INFLUENCE OF TWIST 



As noted in several previous studies ||7|,p[pX|], twist is an important degree of freedom brought by the extension to 
3D. It is well-known from classic studies of elasticity p6| (and everyday experience) that straight rods and ribbons can 
be destabilized by twisting them beyond a certain level. A somewhat similar instability was reported in ref. px| ] in 
numerical simulations of excitable filaments. Beyond a threshold twist, the rotation center line of an initially straight 
twisted filament was observed to adopt an helical configuration. Observations of a similar 'sproing' pO[ ] instability have 
since been made in the related context of the complex Ginzburg-Landau equation vortex lines p7[ . The characteristics 
of the excitable filament sproing instability have however remained somewhat unclear. In the dynamical simulations 
of ref. |]lO|, a single filament turn was imposed in a simulation box with periodic conditions, and the box height was 
varied. A complicating feature of this procedure is that both twist and the available wavelength range are changed 
at the same time. On the theoretical side, the instability fails to be captured by small twist approaches ||l^,|l^ since 
Ip^ the motion of the rotation center is not influenced by twist in this limit (see appendix |d|) . 



^An explanation can be provided by the proximity of the U2 = ijJi point on the meander threshold line (as noted in a particular 
limit in The argument is that i) the meander modes are close to the translation modes when CJ2 is close to lo\ ii) the 

ellipse should reduce to a circle for translations since the translated circular core is circular. This can be explicitly seen from 
Eq. (psh. Namely, Eq. ( p^ gives &x' + iSy' = cos{kzz)[A{a + if5) exp(ia;2i) -I- A* {a* -f i/3*) exp(— iaj2i)]- The explicit inversion 
of Eg. (p6|, p?! ) shows that a + ifi \s, proportional to {vtUi — UtV\) and therefore vanishes when when tends toward the 

translation eigenmode {ut,Vt). In this limit, it only remains the term proportional to a* + ij3* and 5x' + iSy' follows a circle. 
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(a) (b) 




(c) (d) 




FIG. 9. (a): view of the restabilized state in the weakly non linear regime. The parameters are : 
a — 0.684, b — 0.01 and e — 0.025, the size of the simulation box is (128 x 128) x 130 with dx — 0.2 and periodic boundary 
conditions are used. The dashed lines are the trajectories of the tip of the spiral in four regularly spaced horizontal planes. 
The bold line is the instantaneous helical filament. The dash-dotted line is the axis around which the axis of the helical 
filament rotates, (c): black bold circles, projections of the instantaneous filament on an horizontal plane at different times. 
The bold radius in each circle shows the instantaneous filament point at height z = 0. This point trajectory is also shown for 
several periods of rotation (bold dashed-dotted line: evolution between circles 1 and 7; thin gray line: evolution for some time 
afterwards). The pulsation of the axis of the filament is equal to wi = 1.744, the pulsation of the meander in the rotating 
frame is equal to iJ2 — —2.159. These values are to be compared with results of the linear stability analysis : ui = 1.766 and 
iJ2 = ±2.194. (b): view of the restabilized state with the same parameters and same simulation box using no flux boundary 
conditions. The instantaneous fllament has a zig-zag shape and the amplitude of meander varies with z. (d) : black bold lines, 
projections in an horizontal plane of the instantaneous fllament at different times. The trajectory of the spiral tip in a plane 
where the amplitude of meander is maximal is also shown for several periods of rotation (bold dashed-dotted line: evolution 
between fllaments 1 and 7; thin gray line: evolution for some time afterwards). Its pulsation is equal to lji = 1.745 whereas 
the pulsation of the instantaneous fllament in the rotating frame is equal to UJ2 ~ —2.159. The mean distance between the 
instantaneous fllament points and the central (dash-dotted) axis is f?=0.4930 in the (a) (periodic boundary condition) figure 
and R — 0.4944 in the (b) no-fiux case. The core radius of the corresponding 2D spiral is Rq = 0.4833. 



The present approach permits to relieve some of these problems since the twist can be varied from zero to 
large values and a whole range of wavelengths can be examined in the linear stability computations (fc^ is simply a 
parameter which can be given any chosen value independently of r^,). 

We restrict ourselves here mostly to parameter values for which a straight untwisted filament is stable, that is on 
the large a side of the (3D) meander instability region. Fig. |l^ shows the frequency and tip radius for a family of 
twisted scroll wave obtained by increasing from to Tu, at one such parameter point (a = 0.8, b = 0.01, e = 0.025). 
The frequency u!i{tw) increases quadratically at small twist and almost linearly for larger twist values. The quadratic 
behavior at small Tu, is simply obtained by applying first order perturbation theory to Eq. (^, ^) which gives 



= coiir^ = 0) - '"'ft^^^t^^. , + O(r^) (29) 

The direct computation of the matrix element ratio on the r.h.s. of Eq. ( |2^ ) is in good agreement with a direct fit 
of the wi(ru,) curve of Fig. ^ (see caption). Analytic descriptions of the uji{t-u]) curve for larger twist values have 
recently been obtained in the free boundary limit (e 0) both for small core and large core scroll waves . 
The determination of a family of increasingly twisted steady scroll waves permits to determine the evolution of the 
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FIG. 10. (a) : square of the amplitude of the meandering restabilized state as a function of the parameter a. The parameters 
used are b = 0.01 and e = 0.025. The simulation box is (128 x 128) x 130, with dx — 0.2 and periodic boundary conditions 
are used, (b) : comparison of the square of the amphtude of the meandering restabihzed state multiphed by 0.035 and the 
growth rate of the meandering mode Re{a). Same parameter regime as in (a) and a is fixed to 0.68. 




FIG. 11. Frequency uji of the scroll wave as a function of the twist for a = 0.8, 6 — 0.01, and e — 0.025. The small Tu, 
behavior of the pulsation can be well approximated by a'i(ru,) = u}i{t-w = 0) + 0.7205r^ for low values of r^n- For higher values 
of Tu,, a linear behavior of toi as a function of r^, is observed. The first-order perturbation result coefficient ( ]29| ) is 0.7203 using 
the numerically determined rotation eigenmodes of C and of its adjoint. 
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stability spectrum with r„,. The results of such a computation are shown in Fig. |ij for scroll waves of Fig. 111]. As twist 
increase s, the deformations of the translation bands are particularly important. As expected from general arguments 
(section [IB 2), for = 0.2 (Fig. pl]b), the translation bands cr*. 1,0-4.2 are no longer even and related by complex 
conjugation. It only remains the lower symmetry a^-^^^kz) — (Tt,2{—kz)- One can note also that the translation modes 
with R e\(j(k z)\ = stand at kz = iru, and no longer at kz = 0, in agreement with the analytic expression given in 
section |ll Q . When twist is further increased to 0.33 (not shown) a second maximum of Re[a{kzy\ appears near 
kz = 0. The value of Re[a{kz)] is negative at first at this secondary maximum. However, it increases with r^, and it is 
slightly positive at = 0.35 (Fig. |ll|c). The twisted scroll waves are then unstable for a finite range of wavevectors 
near kz = 0. Increasing twist further, enlarges the range of unstable wavelengths and the instability growth rate, as 

shown on Fig. Hi. 

Dynamical simulations reported in section |IV B| show that this twist-induced instability of the translation bands 
correspond to the "sproing" instability of ref. ||10||. Before describing these results, it is worth explaining why the 
instability does not appear around the translation modes at kz = ir^, but a finite wavevector away from them. This 
is a direct consequence of 3D rotational invariance: a twisted scroll wave the axis of which is tilted has the same 
frequency as one the axis of which is vertical. So, a small tilt perturbation should not change the translation mode 
eigenvalues ztiwi to linear order. Therefore, they remain local extrema on the translation bands, as it is observed in 
Fig. (|l^). A direct mathematical proof (based on the same reasoning) is offered in the next section. 
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FIG. 12. Real parts of the rotation band (thin solid line) and the two translation bands (bold dashed and dash dotted lines) 
as a function of the wave number for the same parameter values a = 0.8, b = 0.01 and e = 0.025 and diflerent values of 
twist : (a)r,i, = 0.,(b) = 0.2, (c) = 0.35 and (d) = 0.45. The translation bands have maxima at kz — with a zero 
growth rate. A secondary maximum appears on the translation bands as the twist increases. At a threshold value of the twist 
it becomes unstable at a non-zero value of k^. 



Finally, we find it interesting to show in Fig. ^the twist influence on the spectrum in the "negative line tension" 
parameter regime of Fig. |^, although the untwisted scroll wave is already unstable in this case. Twist modifies the 
spectrum in a way that is rather different from that seen in Fig. ^ It mainly amplifies the instability of the large kz 
part of the spectrum . 
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FIG. 13. Real parts of the rotation band (thin line) and the two translation bands (dashed and dashed-dotted) as a function 
of the wave number for a = 0.44, b = 0.01 and e = 0.025 and different values of twist ; (a)ri„ = 0., (b) — 0.1, (c) = 0.14 
and (d) = 0.19 . The translation bands have a minimum zero growth rate for kz = ±r„. The maximum growth rate of 
the translation bands is increased by twist. The meander modes also become less stable as twist is increased. The change in 
the most unstable band for close to 0.14 appears to be due to hybridization between the translation and meander bands 
like the one observed in Fig. ^ The results of direct numerical simulations show that twist does not qualitatively modify the 
development of the zero twist instability in this regime: a "negative line tension" growth of the filament is observed. 
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A. Helical destabilization and 3D rotational invariance 



In order to demonstrate that the translation modes eigenvalues at kz — ±ru, remain extrema on the translation 
bands, we show that 



da 

dk 



= (30) 



We proceed in two steps. First, perturbation theory is used to compute the eigenvalues of modes close to the translation 
modes on the translation bands. For definiteness, we consider modes at kz — —Tw + 6k, close to the translation mode 
{ut,vt) at kz = —Tyj with a = iuji. Eq. ||) can be written without approximation, 

o{kz)[l\^ ={-i5k^ + 25kT^ + 2iT^5kd^)l^''^^+C^^=.r^ (^^J^ (31) 

Seeking in perturbation, ui — ut + Sui, vi ^ vt + 6vi, one obtains to first order in Sk, 

5a ) + III ) = ^'^^"-(^ + ( ) + ^'^=-'- ( ) (32) 

Multiplying by the left eigenvector {ut,vt) oi Ck.=-T^ for the eigenvalue iuoi and taking the scalar product gives. 



2r,„ 

{ut,{l + id^)ut) (33) 



da 

Thus, the translation modes remain extrema on the translation bands, if 

(ut,il + idM = 0- (34) 

Eq. ( ^^ is a consequence of 3D rotational invariance as we proceed to show. The perturbation corresponding to 
inclining the scroll axis can be found by expressing the inclined scroll in the vertical scroll referential, similarly to 
what was done to determine translation modes (Eq. ( pT| p^). 
One obtains. 



j^f ) = exp[i(cjit + T^z)] f \ = cxp[i{ujit + T^z)] zexp(i0)(9r + -d^) + r^r exp{i4>)d^ f "° 1 



(35) 



and the complex conjugate mode. One can directly check that f j^^), obey the linearized time dependent equations. 
Namely, 

{dt + iT^dlz - dlzV'-^c = i^id^ + rldl^ + "^Id)^! + [dufiu^, vo)u^^„i + 9./(uo, vo)vgl]/e (36) 
dAil = ujid4,v'ill + du9{uo, wo)w™c + d^g{ua, vq)v';^^ (37) 

Eq.|^ shows that (uj„c,Wmc) (Eq. (||))obey 

[/:..=-.„--i]f _)= ) (38) 



The inhomogeneous r.h.s. of Eq. (^8|) comes from the fact that the dz derivative terms in Eq. (^7|) act both on the 
exponential prefactor and on the intrinsic z dependence of Ui„c (Eq. (|35|), while in effect only their action on the 
exponential term is taken into account in Eq. ( |38| ) (through the kz dependence of Ck^ ) ■ 

Eq. (|38| ) directly gives the sought orthogonality relation ( pil ) after multiplying both of its sides by the left eigenvector 
[utyUt) of Ck^=-T„ with eigenvalue itui and taking the scalar product. 
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B. The sproing bifurcation 



In order to study the nonlinear development of the twist-induced instability shown in the parameter regime of 
Fig. [l^c, d, we performed direct numerical simulations of twisted scroll waves using periodic boundary conditions at 
the top and bottom of the simulation box. 

Two kinds of initial conditions were used. The simplest one consisted of two dimensional spirals stacked along the 
vertical (z) axis. The twist was introduced by rotating them around this vertical axis. The main disadvantage of this 
type of initial conditions was that they usually are far from a stationary twisted scroll wave and from a restabilized 
wave when it existed. As a result, reaching the asymptotic attracting state could be very costly in computational 
time. In order to avoid this problem we mainly used initial conditions constructed using results of previous direct 
simulations on a grid with the same values of the parameters and the same horizontal size but of a different vertical 
extension, interpolating linearly the values of the u and v field on the new grid. 

For definiteness, we describe the result for the parameters of Fig. O. 



We first focus on the case when a single turn of twist is initially imposed as in ref. |10|. This case is special for the 
following reasons. On one hand, the previous linear stability results (Fig. |l2|) show that all the potentially unstable 
modes correspond to \kz\ < Tw On the other hand, in a box of height H, the only possible kz values compatible with 
the imposed top and bottom periodic boundary conditions are multiple of 27: / H. Therefore, for a single turn of twist 
= 27r/iJ, there is a single potentially unstable mode in the simulation box and it stands at kz = 0. 

A series of dynamical simulations were performed in boxes of varying heights H . The initial twist was correspond- 
ingly varied from = 0.3 to = 0.5 . 

For low values of twist, the twisted initial state simply evolves toward a straight twisted scroll wave. The instan- 
taneous filament has an helical shape that rotates at a fixed time independent frequency around its vertical axis. In 
each horizontal plane, the wave is seen as a spiral steadily rotating around the helix axis. 

Beyond a threshold twist Tc, the twisted initial state evolves toward a more complex state. The threshold twist is 
estimated to be Tc — 0.352 from the direct numerical simulations. It closely corresponds to the value 0.350 obtained 
from the linear stability analysis for the instability of the fc^ = modes (the instability threshold twist at fc^ 7^ is 
r ~ .345). The asymptotic scroll wave state is shown in Fig. |lj for = 0.381 and it is qualitatively similar for other 
values > Tc- The instantaneous filament takes an helical shape at each time. The axis of this instantaneous helical 
shape is independent of time but its other characteristics vary with time. The point of the instantaneous filament in 



a given horizontal plane (i.e. the spiral wave tip in that plane) closely follows an epicycloidal motion (Fig. 14) and 
the instantaneous filament global evolution can be accurately parameterized as 

X = Ri COs{LJst + T^,Z + 0) + RcOs{u},nt + T^z) 

y ^ RlSU\{uJst + T^Z + (j)) + Rsa\{ujrnt + T^^z) ^ ' 

The pulsation and the radius R\ are close to the pulsation and radius of the stationary straight twisted scroll. 
The pulsation is found to be small compared to (Table. |l|). The radius R is zero at the bifurcation threshold. 
It increases and become comparable to R\ as increases past Tc (Fig. |l5|). As RiOm remains small compared to 
RiUJs, the movement of the spiral tip in an horizontal plane can be described as a rapid rotation movement around 
a slowly moving "mean" point (Fig. ^4|). The bifurcation can thus be described as in ref. ||l^ as a transition in the 
shape of these slowly moving points, the "mean filament" , from a straight shape to an helix of radius R (with the 
same axis as the instantaneous filament). 

The amplitude of the sproing bifurcation can be measured by the radius R of the helical mean filament^. The 
numerical simulations results (Fig. |5|) show that R behaves as ^/tw — Tc confirming the normal Hopf type of the 
bifurcation. The motion of a filament point in an horizontal plane is quasiperiodic with two frequencies (see Eq. |39| 
and Fig. ujm and w^. As reported in Table ||, the frequency ujrn closely agrees with the difference uji — uit.k,=o 
between the stationary twisted scroll pulsation (cji) and the imaginary part of the unstable mode at kz — 0, as 
expected from a Hopf bifurcation in a rotating frame. The other frequency lUs is equal to the twisted scroll pulsation 
u)i at the bifurcation point but departs from it as one moves away from it. 

Some insight into the sproing bifurcation and the value of ujs can be gained by computing the local twist of the 
restabilized scroll wave. The mean filament can be taken as the central curve of a ribbon, one edge of which is the 



instantaneous filament of length. The local twist of this ribbon is equal to (using Eq. (|C2|) and ( C4 ) of appendix O 



can be easily computed as i? = {Rmax + Rmin)/2 where Rmax (Rmin) is the maximum (minimum) distance of the 
instantaneous filament from its axis in an horizontal plane {Rmin should be counted negative value if the spiral tip trajectory 
enlaces the axis of the helices during one rotation period). 
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ioi(Tu.) - Im{at{k^ = 0)) 






r 


LJi(r) 


0.50 


0.053 


0.049 


1.973 


1.91 


0.353 


1.892 


0.45 


0.061 


0.059 


1.943 


1.88 


0.347 


1.889 


0.40 


0.066 


0.066 


1.916 


1.88 


0.347 


1.889 


0.355 


0.063 


0.063 


1.891 


1.88 


0.349 


1.890 



TABLE II. For different values of the imposed twist Tw, values of the pulsation of the rcstabilizcd mean filament ujm, of 
the difference between the pulsation of the steady scroll uj\{tw) and the imaginary part of the unstable translation mode for 
kz = 0, the pulsation of the steady scroll uj\{tw) , the pulsation of the restabilized spiral around the mean filament ujs the local 
twist r of the restabilized state and the calculated corresponding pulsation lo\(t). One notes that los is close to uoi{r) 



L 1 + {Rt^Y ^ ' 

where L is the mean filament length and = 2it/H is the twist imposed on the initial straight scroll wave. Eq. 40| 
shows that sproing decreases the twist. Moreover, when the initial twist is increased the average helix radius R also 
increases so as to maintain the local twist r approximately constant, as shown in Fig. [l5|b. The frequency LOg appears 
to remain close to the frequency of twisted straight scroll wave with this value of the local twist, as shown in Table ||. 

It is interesting to compare the above results with what happens when the initial condition contain n turns of twist 
since the modes with = j2Tr/H, j = I, ■ ■ ■ , n— 1 obey both kz < and the periodic boundary conditions. Analyzing 
moderate values of the imposed twist — n2Tr/ H requires of course to extend the box height H proportionally to n 
and restricted us to n < 5. 

The simplest interesting case occurs when a single unstable mode with k^ ^ Q can develop in the simulation box. 
The linear stability results shows that this can only happen close to the instability threshold when the instability 
growth rate is very low, otherwise the fc^ = mode is also unstable. This is achieved, for instance, for a ~ 0.8, b = 
0.01 and e = 0.025 and n — b initial turns of twist in a box of height H = 613 x dx with dx = 0.15. The single 
unstable mode kz = 27r/ff — 0.069 corresponds to a wavelength equal to the box height H . With this parameter 
choice, a direct simulation shows that the instability develops. The asymptotic state is similar to the previously 
described one for a single unstable mode at kz = 0. The movement of the corresponding instantaneous filament can 
be parameterized using polar coordinates in each horizontal plane by: 

x + iy^ R^e''^^t+^^'- + ii2e+"^2t+{rn,-k.2)z+i, ^ (4^) 

where i?i ~ 0.31 and wi = 1.85 are close to the radius and pulsation of the straight twisted scroll and where R2 ~ 0.1 
and UJ2 = —0.047 are small compared to i?i and uji. The wavenumber kz2 = 0.069 is equal to the single unstable 
wavenumber that can develop in the simulation box. The computation of 102 shows that it is close of the difference 
wi(tu,) — Im{at{kz2)- In contrast to the single turn of twist case, the instantaneous filament shape is slightly different 
from an helix since its radius varies along the vertical axis. The motion can nonetheless be interpreted in the same 
manner by considering that the scroll rotates uniformly around a slowly moving helical mean filament of radius R2 
and pitch 2-1: /{t^ — kz2)- 

The case where several unstable modes of the translation band can develop in the simulation box, can only be 
studied in a box where several turns of twist arc initially imposed. The asymptotic state reached in such a case after 
the instability development, starting from an initially twisted straight scroll wave is shown in Fig. |l6|. It is significantly 
more complicated than in the single unstable mode case and the instantaneous filament shape is rather different from 
a simple helix. 

In each horizontal plane, the spiral tip rotates uniformly around a slowly moving point. The positions of these 
slowly moving points can be numerically computed |^ and used to construct the position of a mean filament. The 
motion of the mean filament is found to be well parameterized in each horizontal plane by 



^This can be done by considering the mean filament as the mean position of the instantaneous filament over a rotation period, 
or by considering it as the instantaneous center of rotation of the instantaneous filament in a plane or by removing the high 
frequency peak in the Fourier spectrum of the instantaneous filament motion. These three methods give very similar results. 
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(a) (b) 

0.8^ 



0.6 




FIG. 14. (a): the bold solid line represents the helical instantaneous filament, the dotted line the mean filament and the 
thin solid lines the quasi-circular trajectories of the instantaneous filament in horizontal planes of equally spaced z. Parameters 

values arc a = 0.8, h — 0.01 and e = 0.025, the simulation box is (128 x 128) x 110 with dx — 0.15, corresponding to Tw = 0.381. 
(b): Modulus of the Fourier transform of the spiral tip complex position {x + iy) in an horizontal plane for the same parameter 
regime. The peak at a; « —0.064 corresponds to the slow movement of the mean filament in the plane while the peak at 
UJ = 1.884 corresponds to the rapid rotating motion of the spiral. The Fourier transform was performed using 2048 points with 
a time spacing of 5t = 0.1969. (c): the thin solid line is the trajectory of the spiral tip in an horizontal plane and the bold 
dashed line is the trajectory of the mean filament in that plane. The mean filament rotates clockwise while the fast rotation 
of the spiral tip is counterclockwise. At a given time, the projection of the mean filament and instantaneous filaments on the 
horizontal plane are circles centered on the helices axis position marked by a star. 



(a) (b) 




FIG. 15. (a)(o) Radius (R) of the helical mean filament as a function of the twist Tw, for a stationary straight twisted scroll 

wave and linear interpolation of for small values of t^, — Tc (thin continuous line) where the computed threshold twist is 
Tc — 0.352. The parameters are a = 0.8 b = 0.01 e = 0.025. (b): local twist of the restabilized filament as a function of Tw 
The dashed line is the line of equation t = Tc, the continuous line is the line of equation t = Tw 
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(a) 



(b) 






FIG. 16. o = 0.8, 6 = 0.1, e = 0.025, the simulation box size is (128 x 128) x 349), the space step is dx = 0.2 and five turns 

of twist arc imposed, which corresponds to r„ = 0.450 (a): The bold soUd hne represents the restabihzed mean filament. The 
instantaneous filament and the quasi circular trajectories of instantaneous filaments in horizontal planes are not shown, (b): 
The dashed and solid bold lines represent a top view of the instantaneous filament at two given times. The three thin solid 
lines are typical trajectories of the spiral tip in horizontal planes. The motion of the spiral tip in each plane is the composition 
of a fast counterclockwise and a slow clockwise rotations, (c): Modulus of the space Fourier transform of the complex mean 
filament position (x{z) + iy{z)). The amplitude of each mode is constant in time while its phase grows linearly in time with a 
constant pulsation a; (fez) shown in (d). For the Fourier modes with an amplitude significantly different from zero, a; (fez) is a 
linear function of kz ■ 
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y R.e^(k'^J''^+'^,t+<t',) (42) 



^y= 



3 



where z denotes the vertical position of the plane and kz^^ are the wavenumbers of the modes that have developed in 



the simulation box. All the observed fci"''' have the same sign as r^. As seen from Eq. 42, the corresponding modes 
have an amplitude Rj which is constant in time and a phase ojjt + (pj which evolves linearly in time. Moreover, 

the pulsations LOj are linearly related to the ojj = uj + kz^^ c, as shown in Fig. |l^. Thus, the mean filament 
parameterization can be rewritten as 



x + iy = e"^*F{z ~ ct) with F{z) = V i^^-e'^'^'"^ (43) 



This explicitly shows that the mean filament deformation propagates as a nonlinear wave of constant shape in the 
vertical direction. It was indeed directly checked that the mean filament shape did not noticeably change at long 
times in our simulation (i.e. for time intervals as long as At = 2000). A direct computation of the twist shows that 
it is almost uniform and that its mean value is significantly lower than the one obtained when only kz — mode can 
grow. In the case depicted in Fig. the mean local twist of the restabilized state is 0.327 whereas in the restabilized 
state when only the kz — mode can be destabilized, the mean local twist is equal to 0.347. 

The stability of the simple restabilized helices was tested in the parameter regime were the more complex state of 
Fig. [l^ existed. To this end, a simulation was first performed with a single turn of twist in a box height chosen such 
that the initial twist was equal to r^, = 0.45. It produced, as already described, a restabilized helix similar to the one 
shown in Fig. |l^. Five copies of this restabilized helix were then vertically stacked. The resulting five turn helix was 
used as the initial condition for a simulation analogous to the one shown in Fig. |l^. It was observed to evolve toward 
a complex state identical to the one shown in Fig. This clearly shows that the simple helices are unstable and 
strongly suggest that the complex state of Fig. ^ is the unique attractor in this parameter regime. We repeated this 
computation with an initial twist of 0.357 and again a similar complex state was produced. 

Systematically varying the initial twist would permit to see whether these complex states directly appear at the 
instability threshold or arise from a secondary bifurcation of stable simple helical states. More generally, further 
studying the filament shapes in longer boxes with more unstable modes would be quite instructive. Both tasks require 
more computer time and power than presently available to us and should be left for other studies. 



V. CONCLUSION 



We have searched to gain and present a general view of scroll wave linear instabilities and of their nonlinear 
developments for a simple model of an isotropic excitable medium. Different types of instabilities have been shown 
to arise depending upon the band of modes to which they belong. These different instabilities have been found to 
develop along different ways and to give rise to distinctly different restabilized states that we have endeavored to 
characterize. 

The negative line tension type of instability has been found to occur in the weakly excitable part of the phase 
diagram and to be strictly linked to 2D spiral drift in an external field. In this respect, it seems worth to try and 
better analyze the mechanisms of spiral drift change. The meandering instability is present in 3D in a larger parameter 
region than in 2D. On the "small core" side of the phase diagram, scroll wave meander in a regime where spirals are 
steadily rotating, as previously noted p^ ]. The long wavelength deformation of the meander band is however not 
directly related to spiral drift. The introduction of twist has been found to induce a deformation of the translation 
bands. Above a threshold twist, this gives rise to the sproing instability which takes place a finite wavevector away 
from the translation mode. This has been shown to be a general consequence of 3D rotational invariance and it explains 
that the sproing bifurcation is not captured by small twist approaches. The bifurcated state arising from the growth 
of a single unstable mode has been found to take a simple helical shape as previously described However, a more 
complex filament deformation has been observed to result from the growth of several unstable modes. Simulations in 
larger boxes appear needed to better analyze these states. 

The present study appears worth extending along several other lines. It certainly is important to investigate how 
the present results extend to more realistic models of excitable media, specially in the cardiac physiology context. It 
will also be quite interesting to see how rotating anisotropy . or spatially varying properties induce the different 

instabilities or interact with them. Deeper insights into the behavior of these complex waves in various situations 
may be gained by developing and analyzing of reduced models reproducing the basic phenomenology here uncovered. 
We hope to be able to report progress in this direction soon. Advances in 3D visualization appear to render 
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possible detailed experimental characterization of scroll wave instabilities and dynamics. We hope that the present 
results will also be seen as a further motivation to undertaking this challenging task. 



APPENDIX A: NUMERICAL METHODS 



1. Determination of the steady state 

Determining the stationary scroll waves consists in solving the nonlinear eigenvalue problem To this end, the 
equations are first discretized on a polar grid of size (iV^ x A^^). This provides a set of 2{Nr x N^) equations (the 
values 

of the r.h.s of (||) on the grid points) with 2{Nr x A^^) + 1 unknowns (the values of the field mq, vq on the points 
of the grid and the pulsation uj). This indetermination comes from the rotational invariance of the problem and can 
be taken care of by setting the value of Uo{Nr, N^) to 0.5. One thus has to find the zeros of 2{Nr x N^) non linear 
equations of 2(AV x iV^) variables. 

This can be done accurately using Newton's method. The linear operator involved in Newton's method is inverted 
using an iterative technique (biconjugate gradient ]35| ) and the starting point is either the result of a direct numerical 
simulation interpolated on the polar grid or the result of a previous computation with slightly different parameters. 
We found that the Newton's method always converged and that the convergence was exponential and allowed to reach 
accuracies of order 10^^ in L2 normQ in a few steps (about 10 which take about an hour using a DEC alpha PWS 
500 Workstation for a grid of 80 x 160 points) whereas direct numerical simulations allowed only an accuracy of order 
unity in L2 norm. 



2. Computation of the linear stability spectrum 

The eigenvectors of Ck^ are denoted by 61,62, - • • and indexed according to the real part of the corresponding 
eigenvalues (Ti > a2 > ■ ■ ■■ In order to accurately compute the eigenvalues of of Ck^ of largest real parts, we have 
used an iterative method proposed and analyzed in pc| ]. It is briefly described in this appendix and some details of 
our implementation are provided. 

The idea of the algorithm is to diagonalize a projection of Ck, in a subspace (approximately) spanned by the 
eigenvectors 61, • • • , 6m corresponding to a given number m of eigenvalues of largest real parts. 

The algorithm proceeds in three main steps. 

The first step consists in creating an appropriate vector xi for generating the diagonalization subspace. A suitable 
choice is to take xi — exp{Ck^t)xo for a generic vector xq. This suppresses the components of xi on the eigenvectors 
on high index [note that these eigenvectors correspond to eigenvalues of large modulus] . The end consequence is that 
truncation at level m produces an error of order exp[t(i?e((Tm+i — cj] on the representation of the i-th eigenvector 
(for i < m) [20I] . In practice, multiplication by the exponential is approximately achieved by computing by iterations 
xi — {I + dtCk^Y"^'^*'XQ, for an arbitrary vector xq, a sufficiently large integer tQ/dt and a sufficiently small dt to 
prevent the time stepping scheme from being unstable (this means that dt is lower than l/max(||Ai||), here dt ~ 10^^). 

The second step consists in generating an appropriate subspace E for diagonalization. This is taken to be the 
space E spanned by (G(£fc^))"xi, n = 0...m — 1 where G is an polynomial, the choice of which is discussed below. 
Explicitly, the computation of an orthonormal base of E and of the matrix A = aij of the projection of (G{C))m on 
E proceeds recursively as follows. Let xi, Xn be the n first element of the base of E and G{Ck^)xn — y„. Then, 
at the next step we compute: 

„ _ Vn " '^i=l...rXyn-Xi)Xi (\-\\ 
Xn+1 — II ^ / \ II ■ l-^-'-J 

The construction ends at step to. This building of an orthonormalized base of E presents the advantage of decreasing 
the contribution of the first eigenmodes in the elements of higher order of the base. Otherwise, this contribution would 
be dominant and would result in a lower accuracy in the computation. In our case, the scalar product is defined by 
a discretized version of the standard scalar product (p^): 



'We define ||«, w||2 = + ''^Xj) 
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) = ^ ^ Xi{ir, ie)yj{'ir, ie) ir drd9dr 



(A2) 



The third step consists in the diagonahzation of an appropriate truncation of {G{Ck,))- Since E is not invariant 
under the apphcation of {G{Ck^)), the orthogonal projection Gm of G(Ck^) onto E is considered. Its matrix elements 
are: 

9i,j ^ {xi,yj), (i, j) e {l...r7i} X {l...m} (A3) 

The obtained matrix is diagonalized and both its eigenvectors and its eigenvalues are computed. It is checked at the 
end that the leading eigenvectors of Gm are good approximations of the leading eigenvectors of £k^ ■ 

The choice of the polynomial G is of some importance. Indeed, the simplest and computationally most efficient 
choice would be to take G{X) = X. This would result in the amplification of the contribution of the eigenmodes of 
of large index [i.e corresponding to eigenvalues of large modulus] and, Ck^ being ill-conditioned, would prevent 
the success of the method. We have used G{X) = (1 + dtxy^/''-*, with ti chosen large enough to make G{Ck^) 
differ significantly from the identity (a typical value was ti = 0.5). Despite the increase in computational cost, this 
choice significantly increases the accuracy of the method by decreasing the contribution of the eigenmodes of negative 
eigenvalues. Using m = 50, the 10 most unstable eigenmodes and eigenvalues are obtained with a good accuracy: 
\\{a~CkJ{ui, vi)\\2 < 10-6- 

3. Direct numerical simulations 

Direct numerical simulations of three dimensional excitable media were performed using a forward Euler explicit 
timestepping scheme . The diffusion operator was evaluated using finite differences and a 19 points formula js^, 

+ 2(Mi+lj_fc + Mi-lj\fe + Wij-+l,fc + Mij-l,fc + + Uij^k-l) 

+ Wi+lj + l,fe + "fij+l j-l,fc + Wi+lj,fc+l + ■"i+lj,fc-l 
+ Ui-l,i+l,fc + Ui-l,j-l,k + + "Uj-lj.fc-l 

+ + + Uij + i^k~l + + Ui,j~l,k-1 (A4) 

This method has two main advantages compared with the classical 7 points formula. First its stability limit allows 
greater time steps and therefore the computing time in spite of the additional operations involved in the Laplacian 
evaluation is found to be 1.3 times smaller. Second, the error made when evaluating the diffusion operator is isotropic 
at the dominant order(order dx'^), whereas with the 7 points formula it depends of the grid orientation. 

No-flux boundary conditions n.Vw = are imposed on the vertical sides of the box and either no-flux or periodic 
boundary conditions are chosen on the top and bottom boundaries. 

The position of the instantaneous filament in the horizontal planes is computed as the intersection of a u = 0.5 and 
a, V — 0.75 (0.5a — b) isosurfaces. 



APPENDIX B: SPIRAL DRIFT IN AN EXTERNAL FIELD 



In the presence of a small external field E, the spiral rotation center drifts at a constant velocity proportional to 
the field magnitude but at an angle with the field direction pTj-pOj as given by Eq. (|l8|) of the main text. The drift 
coefficients a|| and a± have been computed in the fre e bo undary limit (e 1) both for small core and large core 
spirals ||lj]. We derive here a general formula (Eq. ( |B10| ) below) for the drift coefficient a = a\\ + ia± as a matrix 
element between the translation eigenvector {ut,vt) of C and the corresponding right eigenvector [ut,vt) oi C. 

With an added external field E, Eq. (|l|, ^) read 

dtu = V^u + f{u,v)/e - E.^u, (Bl) 
dtv^g{u,v). (B2) 

We choose a coordinate system with the x-axis along the field direction and corresponding polar coordinates (r, 0). 
We analyze the motion of a counterclockwise rotating spiral which is stationary,for _E = 0, in the rotating referential 
(r, (/)) with (f) = 6~ujit. Since dx — cos{d)dr ~ sm{d) / rde , the field term E.Vu can be written in the rotating referential 
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a 


r.h.s of Eq. (|B10| ) 


a|l — ia± 


0.44 


-2.05 -0.78i 


-1.97-0.84i 


0.62 


3.5 - 0.47i 


3.4 -0.42i 


0.7 


1.59 - O.SOi 


1.62 -0.83i 



TABLE III. For b = 0.01, e = 0.025 and several values of a, the drift coefficients given by Eq. (BIO) computed using the 
eigenmodes of jCk^ and of its adjoint and the drift coefficients a« — ia± measured in direct numerical simulations. 



EdxU — —E[cos{(j) + ujit)drU — sin(0 + LUit)/rd,pu\, 

= —E /2[e-K.p{iujit) exp(j(/))(9rU + id^tu/r) + c.c. 



(B3) 



As will be seen below, secular terms appear when the E term on the r.h.s. of Eq. (Bl)) is treated in perturbation. 
Their origin is, of course, the induced spiral drift. Anticipating this phenomenon, we suppose that the spiral rotates 
steadily in the frame with coordinates {r,6) which drifts with respect to the lab frame with coordinates {x,y). That 
is. 



and 



X ^ xo{t) + r cos(0 + Wit), 
y ^yo{t) +r sin(^ + uit), 



dt\r,(i> = dt\x,y + i^idcf, + xodx + yody. 



(B4) 
(B5) 



(B6) 



The supplementary terms proportional to xq and j/o should be chosen to cancel the unwanted secular terms. This 
determines the spiral drift. Explicitly, we rewrite Eq. (^ ) using the rotating coordinates as 

dt\x,y = dt\r,4, - LOid^ ~ l/2[(io - iyo) exp(iwit) exp(i0)(5r + id^/r) + c.c.]. (B7) 

Substitution of this formula in Eq. (Bl .B^) and linearization under the form u = uo{r, cj)) + Up expiiujit) + c.c, v = 
vo{r, (p) + Vp exp{iuJit) + c.c. gives, 



{^u;^ -C)( y ]- l/2(io - *2/o) 



-E/2 



(B8) 



Multiplying by the left eigenvector C of eigenvalue icui would show the need for secular terms if we had not introduced 
the drift terms. Here, however it simply determines the drift as 



{Ut,Ut) 



io-iyo^E — 

{Ut,Ut) + {Vt,Vt) 

Equivalently, this gives the sought formula for the drift coefficients 

{ut,ut) 



a 



{ut,Ut) + {vt,Vt) 



(B9) 



(BIO) 



APPENDIX C: TWIST AND WRITHE OF A RIBBON 



We recall the definition of quantities associated to closed ribbons and some useful mathematical properties for 
analyzing the "sproing" bifurcation. In particular we give a mathematical definition of the twist and its value in the 
case of a uniformly twisted ribbon with an helical central curve. 

The local twist of a ribbon is classically defined as ( |^ ) 

r=(^(plAp).r (CI) 

where t is the unit tangent vector to the mean curve of the ribbon, s is the curvilinear coordinate along the mean 
curve and p, the unit vector perpendicular to that curve that directs the line intersecting one of the edges of the 
ribbon. The twist measures the spatial rotation rate of the edges of the ribbon around the mean curve. 
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For a closed ribbon, the linking number Lk is the integer which measures the entanglement of the two edges of the 
ribbon. Lk is a topological invariant which is constant under a continuous deformation of the ribbon (as long as it 
does not intersect itself). 

The linking number is related to the integral of twist by a formula |^o[] (which has been popularized by its use in a 
molecular biology context) 

Lk^Wr + ^ f Tds (C2) 



2tt _ 

The writhing number Wr depends only of the mean curve r(s)of the ribbon and is equal to: 

Att J J \\r-r'P 

The tangent vector to r(s) traces a closed curve on the unit sphere as r(s) goes around the ribbon. The writhing 
number is also equal, up to an even integer, to the area enclosed on the unit sphere by this closed curve divided by 

An example of interest is the writhing number of a single turn of an helix of radius R and pitch H (linking the 
two free ends by a non self intersecting planar curve to obtain a closed curve). The tangent vectors curve encloses 
a spherical cap on the unit sphere of normalized area equal to (1 — cosO) where 9 is the angle made by the tangent 
vector with the vertical axis. This gives the writhing number (up to an even integer which is seen to be zero by using 
the continuity of Wr and the fact that the writhing number of a straight line is equal to zero) , 

Wr = l = (C4) 

^1 + {2^R/Hf 



APPENDIX D: LINK WITH AVERAGED EQUATION 

In ref . p6[ , equations were derived for the motion of the mean scroll filament and the evolution of twist for a weakly 
twisted and weakly curved scroll wave. It was subsequently noted fl^ that many coefficients in ref. original 
equations were identically zero and that only four non-trivial ones remained to be determined. This approach has 
recently been extended to take into account fiber rotation anisotropy ||33|] . 

In this appendix, we find it of some interest to explicitly relate the averaged equations of ref. ||lq;|2[ to the 



computations that are performed in the main part of the present paper. One coefficient in the equations of ref. [g6 12 
is given by the quadratic scroll rotation frequency dependence at small twist (Eq. (|29|)). Unsurprisingly, the three other 
ones are given by the curvature of the rotation and translation bands around the corresponding symmetry eigenvalues. 
This explicitly confirms that the sproing instability cannot be captured in the limit considered to derive the averaged 
equations since the instability takes place a finite distance away from the translation symmetry eigenvalues on the 
translation bands. 

We provide a simple-minded derivation of the equations ref. ||l6|,|l2| using a Cartesian frame instead of the more 
sophisticated intrinsic mean filament coordinates used in ref. |ll6| . This limits us to consider a weakly inclined filament 
but we can proceed very similarly to the derivation of spiral drift in appendix ^ and the extension to non-isotropic 
or non-homogeneous medium is straightforward (but not considered here). 

We denote the fixed Cartesian coordinates by (x, y, z) and by (r, 0, z) the cylindrical coordinates of a frame rotating 
at the 2D spiral frequency and centered on the mean filament position (xo(t, z), yo{t, z)), 

X — XQ{t, z) -\- rcos[0 -I- uit + ^{t, z)] (Dl) 
y = t/o(<, z) + r sin[0 -I- uit + ij{t, z)] (D2) 

The corresponding relations between the time and vertical (z) derivatives in the two referentials arc thus, 

dt\r,4, = dt\x,y + ^id^ + dtil> + dtxo + dtyo dy (D3) 
9^ = i9z U,a + d^ij + d^xo + d^yo dy (D4) 

These relations can be rewritten using dx +idy — e-xp[i(f) + iuiit + iilj{z,t)][dr + id^/r] and introducing wq = xq — iyo, 
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dt\x,y = dt\r,^ - Loid^ - 9*^9^ - ^ [dtw^ e'^+^i (9, + id^/r) + c.c] (D5) 
dA^,y = - d^^d^ - i [a.uio e'^+*'^i*+^'^ {dr + ^a^/r) + c.c] (D6) 

+ i [(a,u;o)' e2('"i*+*^) (e'*(9, + id^/r)f + c.c" 
+ \\d.wo\^ylD - \ [9L«^oe**+''^^*+''^ (a, + id^/r) + c.c] 

- + 9.^ [9.u;oe'*+'"^*+^'^ {d% + zcJ^^) + cc] (D7) 

With Eq. (D5, D7) , the governing reaction-diffusion equations (0J|) become 

{dt - ujid^ - VId) u - f{u, v)/e = dti^d^u + i [dtwo e"f'+'''''+'^'{drU + id^u/r) + c.c] + dl,u 

+ - aa.^a^^u - [a,i«o e'*+'"^*+^'^(a2,^. + idl^u/r) + c.c] 

+ i [(9.u;o)' e2('"i*+^'^) (e'*(a, + id^/r)f u + c.c 

+ ^|9.u;oPv2^7/ - i [a2,«;oe^*+'"^*+*'^ (9.^/ + id^u/r) + c.c] 
- d^u + [a,u;oe'*+^-i*+^'^ (a^^u + idl^u/r) + cc] (D8) 
{dt - Luid^) V - g{u, v) = dt^d^v + i [dtwo c'^+''^^'+''' [drv + id^v/r) + c.c] (D9) 



For a wealcly curved and wealcly twisted scroll wave, the r.h.s. of Eq. (D8,D9) can be treated in perturbation starting 
from the 2D spiral fields {uq,vq). At first order, the inhomogeneous r.h.s are a superposition of time independent 
terms and of terms oscillating at frequencies wi and 2wi. One can therefore seek (m, v) in perturbation as 



Wo 
Wo 



(2(^1) 

^1 

(2i^i) 



+ 



^1 

„('^i) 



e'"i*+*'^ + c.c. 



.,(0) 



(DIO) 



(2u;i) (2a;i) , (i^i) ((^i) 



where (u^ , , ("i 
order perturbative corrections. 
The 2a;i-functions obey 



and {uf \ v'f''^) are time-independent fmictions characterizing the three different first 



(2wi) 



'^{dr+id^lr)yu 







(Dll) 



The operator {i2uji — L) is invertible and the 2uJi functions (w^'^^'', dJ^'^^'*) can be determined. They simply describe 
the inclined circular core which is viewed as elliptical in the chosen xy coordinates (these terms are absent in the 
filament coordinates used in ref. where in effect dzWo is zero). 

On the contrary, the ui and constant functions arise from resonant forcing and can only be determined when 
solvability conditions are verified. 

The LOi functions obey, 



(D12) 



Since (ut,vt) is an eigenvector of £ with eigenvalue iuJi, Eq. (D12 ) is solvable only if its r.h.s has no component on 
this eigenvector. More explicitly, one obtains by multiplying Eq. (D12) by the associated left eigenvector. 



dtWQ - dl^wo — 



\Ut,Ut) 



Ut,Ut) + {Vt,Vt) 







{Ut,Ut) 



(D13) 



The last term on the l.h.s should vanish by rotational invariance since simply inclining a twisted scroll cannot induce 
its drift. This explicitly follows from Eq. (pj) since 
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{ut, e'^idl^u^ + idl^u^/r)) = (ut, (1 + id^ = 0- 



(D14) 



Eq. ( pl3D [without the last term] is the equation of motion for the mean filament obtained in ref. ||l6|,|T2|. In 
the limit considered, the mean filament motion is independent of the scroll twist and only depends on the filament 
curvature (d^^wo) with a coefficient which gives both the small dependence of the translation bands (Eq. (|2^)) and 
spiral drift in an external field (Eq. (BIG)). 



The time independent component {i^,v^'') of the first correction (DIO) obey 



C 



,.(0) 







:\dzWo\ 



(D15) 



Again since the rotation mode is an eigenvector with eigenvalue zero of C, this equation can be solved only if 



U^, Ucj,) + (W0, V^) 



(a.V^) 



12 (m0,V2d"O 

2 {uij,,u^i 



= 



(D16) 



Rotational invariance again implies that the last term on the l.h.s. of Eq. (D16) vanishes (since simply inclining a 
scroll wave cannot change its rotation frequency) as explicitly shown below (Eq. (D18). The remaining Eq. (D16) 
is the equation obtained in ref. 0,11^ for the scroll twist dynamics. The coefficient of {dzip)'^ simply describes the 
twist dependence of steady scroll rotation frequency (Eq. (P9|)) while the coefficient of d'^^^p governs the small fc^ 
dependence of the rotation band [Eq. ( p3| ) with the subscript t replaced by (/)]. 

The computatio ns reported in the present paper permit to explicitly evaluate the four real coefficients which appear 
in Eq. ( D13| , D16). For example, for a ~ 0.8, b = 0.01, e = 0.025, the reported results give 



{Ut,Ut) 



{Ut,Ut) 



= 0.8842 - 0.662, 



= 0.578, 



-0.7203. 



(D17) 



More generally, the two coefficients in Eq. (D16) do not change sign as one traverses the different regions of the 
phase diagram: a small twist increases the scroll rotation frequency and a small modulation in the z-direction increases 
the stability of the rotation band modes. The complex coefficient of Eq. (D13) is directly linked to the scroll wave 
line tension stability/instability and its real part change sign as 2D spiral drift. 

Finally, it may be worth comparing the above derivation to that of [ p^ . Since the l.h.s of Eq. { pWj ) is a superposition 
of terms with different time dependences, this is also the case of its solution and many coefficients formally introduced 
m ref. do not even appear in our derivation, as previously noted in [|l^ (in a slightly different formulation). On 
the other side, we have chosen a simple parameterization for which rotational invariance is not manifest, in contrast 
to ref. jl^]. This forces us to explicitly show that coefficients which do not appear in ref. [^ vanish. 

We conclude by showing that indeed, 



{u^,VIjjUq) = 0. 



(D18) 



This is a simple consequence of the transformation property of the reaction-diffusion equations (|[ ^ under dilation 
Namely, {uo{r{l + rj), 0), voir{l + rj), 0) is a solution of (|5[ ^) with V^d replaced by 1/(1 + rj) 
77 <C 1 gives the infinitesimal version of this transformation 



^Vjjj. Expanding for 



rdrVQ J V 



(D19) 



which can also be directly checked by differentiating Eq. (|5|, |6D with respect to r. The multip licati on of Eq. ( D19 ) on 
both sides by (u^, w^), the left eigenvector of C of eigenvalue zero, gives the desired identity (D18). 
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